Privacy-Preserving Framework for Multi-Institutional Medical Time-Series Analysis via Homomorphic Encryption: Design and Development Study

Introduction

Background

The integration of massive, multiparty patient data is crucial for developing high-quality medical AI. Compared to single-site datasets, those aggregated from multiple institutions exhibit greater diversity in patient demographics, clinical practices, and health care conditions. This enables researchers to obtain sufficient samples to train robust models for specific tasks, such as studying rare diseases [1]; achieve results with higher reproducibility and generalizability [2]; accelerate the translation of research findings into clinical practice [3]; and discover novel patterns that are undetectable in any single dataset [4]. As a result, such integration significantly promotes the progress of medical data analysis and makes their applications widespread in many real-world clinical scenarios.

However, the sharing of such data is severely constrained by various privacy regulations, such as the HIPAA (Health Insurance Portability and Accountability Act) Privacy Rule in the United States and the GDPR (General Data Protection Regulation) in the European Union [5,6]. These regulations, coupled with the severe repercussions of privacy breaches—including reputational damage, substantial financial penalties, and health care fraud—make medical institutions highly reluctant to share their sensitive patient data [7]. Consequently, there is an urgent need for secure multiparty medical AI analysis protocols that can navigate the nontrivial trade-off between practical utility and robust security.

Prior Work

Most existing privacy-preserving AI protocols are designed for, or evaluated on, cross-sectional data. This paper, in contrast, addresses the challenge of analyzing longitudinal patient data using recurrent neural networks (RNNs) in a privacy-preserving manner. Longitudinal patient data, which records patient visits, test results, and treatments over time, provides more invaluable insights into the relationship between temporal patterns of features and events of interest [8]. RNNs are particularly well-suited for mining such data due to their ability to model nonlinear relationships and high-dimensional temporal dependencies. They have proven essential in advancing personalized medicine through disease progression modeling [9], outcome prediction [10], and treatment recommendation [11]. Despite their value, longitudinal datasets are often harder to collect and tend to be smaller in scale than their cross-sectional counterparts. This scarcity frequently necessitates pooling data from multiple providers to train robust models, a process thwarted by the distributed and sensitive nature of the data itself.

Nowadays, numerous privacy-preserving deep learning solutions have been proposed for architectures like multilayer perceptrons [12], convolutional networks [13], and generative adversarial networks [14]. However, few have focused on secure RNN-based protocols. Among them, many used federated learning and have been applied in various tasks [15-18], allowing parties to share only local intermediate results like gradients, thus enabling the training of predictive models without direct original data exposure [19]. However, several studies have shown that such aggregated information can still reveal parties’ private data [20,21], leading to privacy compromise risks even with large batch sizes and noisy gradients [22-25], significantly hindering the adoption of these solutions. Furthermore, some recent works have also found that RNN-based algorithms are more vulnerable to attacks targeting aggregated information (eg, membership inference attack) and less well-protected with differential privacy than other common neural networks, making the privacy-utility trade-off even more challenging [20,26,27].

In contrast to federated learning, cryptographic techniques like secure multiparty computation (SMPC) offer stronger cryptographic security guarantees. However, SMPC frameworks often incur substantial communication and computational overhead that escalates with the number of participating parties. This scalability challenge remains a critical bottleneck for large-scale applications. Despite ongoing optimizations, most SMPC-based RNN implementations are still confined to 2-party or few-party settings [28,29], which is often insufficient for multi-institutional medical collaborations.

Homomorphic encryption (HE) enables computation on encrypted data, making it a naturally appealing approach for privacy-preserving analytics. However, its practicality is constrained by a limited set of supported operations—primarily additions and multiplications. This fundamental limitation necessitates the approximation of activation functions using polynomials, introducing a critical trade-off among approximation accuracy, computational efficiency, and noise management. Consequently, nearly all prior attempts at HE-based machine learning and deep learning analysis have incurred more or less accuracy loss. Some try to maximize efficiency but at the cost of a more significant loss in accuracy. For example, Podschwadt and Takabi [30] proposed a noninteractive RNN prediction protocol using parallel RNN blocks, which minimizes the number of communication rounds at the cost of a 5%-9% F₁-value drop; Tosevski and Gulak [31] proposed an RNN quantization procedure to reduce the required HE computation, but their top-5 accuracy of secure speaker identification decreases from 0.9252 to 0.7805. While some other HE-based solutions achieved similar time-series model performance to nonsecure counterparts [32-34], it is crucial to note that a majority of these works focus solely on secure inference—performing predictions on an already-trained model. This approach presupposes the existence of a central, nonencrypted model, which is precisely the scenario that multiparty collaboration aims to avoid. The more challenging and impactful problem of secure training—learning a model from scratch on distributed, encrypted data without ever decrypting it—remains largely unaddressed for RNNs. A limited number of studies have ventured into secure HE-based RNN training, among which Paul et al [35] proposed a collective learning protocol on multiple locally trained RNN models. However, it is essentially a logistic regression training (only fully connected layer parameters change during the protocol) and requires direct sharing of local parameters, which may lead to privacy breaches. Another solution proposed by Sav et al [36] is more similar to ours. They extended ring learning with errors (RLWE)–based HE to the multiparty RNN training scenario by the use of threshold HE, where the decryption key is shared among participants. This introduces a critical new vulnerability: recent attacks have demonstrated that partial decryption shares can leak information about the secret key [37]. The primary countermeasure—injecting significant noise during the partial decryption phase—significantly erodes the precision of the underlying plaintext message, often rendering the model unusable, thus creating an unresolved tension between security and utility in multiparty HE-based training [38]. This unresolved tension between security and precision constitutes a major obstacle for practical multiparty HE-based training.

Study Objectives

To overcome these limitations, this paper proposes a privacy-preserving medical AI analysis framework for distributed longitudinal patient data based on approximated threshold HE [39]. The training process is performed among multiple data providers and a researcher, where all shared intermediate results are encrypted or masked, and only the final model parameters are exposed to certain parties. To achieve practical efficiency without compromising security, we introduce several innovations: (1) a novel matrix multiplication scheme under HE for fewer multiplication depth possessions with very little waste of ciphertext slots; (2) a secure ciphertext refresh strategy that further reduces depth requirements, preserves ciphertext precision, and is proven secure under the indistinguishability under chosen plaintext attack with decryption oracles (IND-CPA^D) security definition against a majority of dishonest passive adversaries [40,41]; and (3) the use of low-degree polynomial and piecewise functions to approximate activations effectively. Additionally, we use a sliced network structure to enhance convergence and accuracy for long sequences [42]. We implement and evaluate our framework on 4 real-world longitudinal medical datasets to assess its security, scalability, and predictive performance relative to nonsecure centralized training.

Methods

Ethical Considerations

This study involved secondary analysis of 4 publicly available, deidentified datasets: one from PhysioNet Challenge 2012, one from PhysioNet Challenge 2019, and two from MIMIC-III (Medical Information Mart for Intensive Care–III).

The MIMIC-III database was approved by the institutional review boards (IRBs) of the Beth Israel Deaconess Medical Center and Massachusetts Institute of Technology [43]. The PhysioNet Challenge 2012 dataset was collected from MIMIC-II, which was approved by the same IRBs of the MIMIC-III [44]. The PhysioNet Challenge 2019 dataset was derived from 3 hospital systems: Beth Israel Deaconess Medical Center, Emory University Hospital, and a third, unidentified one with approval from the appropriate IRBs [45].

This study used only publicly available, deidentified data and did not involve direct contact with participants. Therefore, no additional ethics approval or informed consent was required for this secondary analysis.

Overview of the Proposed Protocol

The architecture of our framework is shown in Figure 1, where two types of entities are mainly involved. One is data providers who hold longitudinal patient data and are willing to share it if no data leakage occurs. The other is a researcher who wishes to obtain a predictive model for various purposes like research or prediction while avoiding any exposure of model parameters to unauthorized parties. We assume that all medical data are transformed into a common format (eg, the OMOP [Observational Medical Outcomes Partnership] common data model [46]) to obtain correct multiparty research data.

As for adversary settings, given the number of involved parties z, up to z − 1 dishonest, passive adversaries can read the secrets of other corrupted parties while keeping the defined protocol running correctly. For the security of approximated HE, considering the revealing of partial decryption in threshold HE is similar to making a decryption query in normal single-key HE, we add Gaussian noises during ciphertext decryption and refresh according to Li’s work to defend IND-CPA^D attack [40], which was recently found to be a threat in both approximated and exact HE schemes, like Cheon-Kim-Kim-Song (CKKS) and Brakerski/Fan-Vercauteren (BFV) [47].

The pipeline of HE-based multiparty RNN model training includes the following steps. First, all parties initialize, encode, and encrypt their shares of model parameters, then secretly sum them into ciphertexts of initial model parameters, which will be broadcast to all data providers. Second, two steps are performed repeatedly until all encrypted local gradients are obtained: (1) each data provider conducts the forward and backward propagation process of the RNN model using its local data and ciphertexts of model parameters till the exhaustion of multiplication depth and (2) all parties together perform the online ciphertext refresh protocol to recover the multiplication depth. Third, all encrypted gradients are summed, masked, and partially decrypted (done by all data providers), and only the researcher gets the masked plaintext gradients. Lastly, the researcher secretly updates its model parameters without knowing the accurate gradients, which will be used in the next iteration. Furthermore, when the training process ends, the encrypted model parameters should be refreshed with a scale-up factor to counteract the impact of partial decryption noises on the accuracy of the decryption. We provide more details about ciphertext matrix multiplication, refresh, decryption, and sliced-RNN training below.

‎

Approximated HE and Optimization Techniques

HE is a type of encryption that allows computations on ciphertexts whose decrypted results can match the results of equivalent operations performed on plaintext. In our proposed framework, to perform secure model training, we use the RLWE-based HE scheme called CKKS that supports addition, multiplication, and vector rotation over the real number [48] (see more details in Multimedia Appendix 1 in the section “Threshold CKKS Scheme”). Since most of the operations over ciphertexts (eg, multiplication and vector rotation) heavily rely on polynomial multiplications, several optimization techniques are applied to maximize the efficiency of the CKKS scheme. These techniques include (1) fast polynomial multiplication based on the number theoretic transform [49], (2) the residue number system (RNS) representation that significantly reduces the computational cost of multiplication between two big integers [50], and (3) the batch encoding technique that makes a single ciphertext contain thousands of values, where any operation can be performed in a single instruction, multiple data (SIMD) manner [51]. Besides, we use a recent error reduction technique by Kim et al [52] so ciphertext rotation can be performed in smaller encryption parameters, and the threshold encryption technique is also integrated for multiparty applications [39]. Under threshold HE, each party holds a secret key, from which all collective keys in public are generated. When computation tasks start, all values will be encrypted by a collective public key, calculated by collective switching keys, partially decrypted by most of the parties, and finally decrypted by the remaining parties so that we can achieve access control of encrypted information and security against passive adversaries with a majority of dishonest parties.

Then, to optimize matrix multiplication over ciphertexts, we propose a novel matrix packing strategy that is suitable for batch and mini-batch gradient descent in common deep learning algorithms, as shown in Figure 2. Suppose an input matrix x and a model parameter matrix β whose sizes are d×v and d×d, respectively (the batch size is usually quite large to make full use of ciphertext slots, meaning that in many cases v>d), we copy the model parameters row by row and encode the input layer as follows. Let m_i denote the ith element of an encoded ciphertext, x_{i_j} denote the element of the input layer in the ith feature and jth sample, l denote the slot count of one ciphertext, and we assume that l/d² is an integer, we have:

${\displaystyle {\displaystyle \begin{array}{c}{x}_{i\mathrm{\_}j}=m_{\frac{l\cdot i}{d}+j}i\in \left[0,d\right],j\in \left[0,\frac{l}{d}-\frac{l}{d^2}\right]\end{array}}}$(1)

Compared to the existing encrypted matrix multiplication protocol proposed by Jiang et al [53], our protocol has the following two advantages: (1) the number of ciphertext rotations is reduced from $3d+5\sqrt{d}$ to $2d+5\sqrt{d}$ with the baby-step giant-step algorithm; (2) the occupied multiplication depth is reduced from 3 to 2, which is much more important, meaning that we can perform more operations with the same encryption parameters or the same operations with smaller encryption parameters. The only drawback of our strategy is reducing the usage rate of ciphertext slots from 1 to $1-d^2/l$, which is acceptable since d is usually small and one ciphertext can contain at least thousands of values. By contrast, Rizomiliotis and Triakosia [54] recently proposed an encrypted matrix multiplication protocol with two multiplication depths but at the expense of a vast waste of ciphertext slots, which is not practical in data-intensive scenarios. Furthermore, when the size of the model parameter matrix becomes quite large, our packing strategy can be combined with the divide-and-conquer optimization scheme to minimize the impact of slot waste. Specifically, we can equally split the matrix of size d² into much smaller matrices of size d'² to increase the usage rate of ciphertext slots from $1-d^2/l\mathrm{t}\mathrm{o}1-d^{\prime 2}/l$.

‎

In the CKKS scheme, one ciphertext’s multiplication operation is limited. Once exhausted, the ciphertext refresh operation should be executed to recover the multiplication depth. Typically, it is done by a time-consuming technique called bootstrapping, whose efficiency is severely dragged down by many linear transformation operations, the fitting of the mod function using high-degree polynomials, and the need for large encryption parameters. Instead, considering the trade-off between computation and communication, we propose an online protocol to securely refresh ciphertexts without plaintext exposure, as shown in Figure 3. Initially, all parties but one generate partial decryption shares and add smudging noises. The noises should be sampled from a relatively large Gaussian distribution to avoid a full key recovery attack (which primarily exploits the approximate decryption results and is implied by IND-CPA^D security) and the leakage of plaintext simultaneously [40]. Next, all shares are sent to the remaining party. Then, according to the received shares, the remaining party fully decrypts the ciphertexts, obtains the masked plaintexts, and transforms the plaintexts as wished. Simultaneously, the same transformations are performed on all generated noises, which will be encrypted and sent to the remaining party later. Finally, after receiving the transformed encrypted noises, the remaining party performs addition among them and masked plaintexts to recover the original information secretly. By allowing transformation on plaintexts, we can perform matrix transpose and the first step of matrix multiplication out of the HE environment, thus eliminating many ciphertext rotations and further reducing the multiplication depth required for matrix multiplication from 2 to 1, which greatly reduces the noise generated during HE computations. Besides, an extra round of key generation should be performed after a certain number of ciphertext refreshes to maintain the IND-CPA^D security, which is not time-consuming if the encryption parameters are small.

‎

Since the exact variance of the smudging noise is contingent upon the noise inherent in the decrypted ciphertexts, a preliminary local simulation is imperative. Prior to the multiparty RNN analysis, each data owner conducts a local training simulation on a small subset of their plaintext data. This simulation strictly adheres to the agreed-upon data preprocessing and model initialization protocols. Its objective is to determine the maximum output bounds for each layer, which in turn informs a static estimation of the noise growth per layer during subsequent homomorphic computations.

Like the online ciphertext refresh protocol, when the model training process ends, the final model parameter decryption also includes noise generation, partial decryption, and full decryption, with the researcher acting as the remaining party. Without adding transformed noises, we refresh model parameter ciphertexts with a large scale-up factor before decryption to avoid the impact of decryption noises on decryption accuracy.

Privacy-Preserving RNN Model Training

Here, we describe how to perform RNN model training securely in two aspects: (1) the low-degree polynomial and piecewise fitting of activation functions and (2) the network structure used in conjunction with the fitting polynomials for better accuracy and convergence without too much computational overhead.

In nearly all HE-based machine learning and deep learning algorithms, a major challenge is the evaluation of nonpolynomial activation functions. All fitting strategies (eg, Taylor series, least squares error, and minimax technique) can only guarantee a reasonable approximation level within limited intervals. Once the input exceeds them, the error will grow dramatically, leading to the possible failure of model training, and things are usually even worse in RNN-based solutions due to the exploding or vanishing gradients. Therefore, selecting appropriate activation functions that are suitable for training tasks and easy to approximate is essential.

For the hidden layer in RNNs, tanh is a commonly used activation function, but its curve is quite steep if the input is close to 0 and very flat elsewhere. When fitting functions with such property, a polynomial with a higher degree should be applied to maintain a sufficient approximation interval at the cost of considerable computational overhead, or the accuracy loss of analysis results is likely to be significant, as demonstrated in previous work [30]. Thus, we choose the rectified linear unit (ReLU) function and use a quadratic function f₂(x) and a 4-degree function f₄(x) to approximate it with an intercept (0, 0) and least squares error, as shown in Figure 4.

${\displaystyle {\displaystyle f_2(x)=0.5x+0.0042x^2}}$(2)

${\displaystyle {\displaystyle \begin{array}{c}{f}_4\left(x\right)=0.5x+0.0109x^2-0.000000647x^4\end{array}}}$(3)

The occupied multiplication depths are only 1 and 2, and the approximation intervals are (−150, 150) and (−100, 100), which are obtained empirically and should be enough for model training on longitudinal patient data with most inputs within them. Even if a few inputs are outside the interval, it will not have a severe impact due to the slow relative error growth of the low-degree function compared to the fitting of sigmoid and tanh (usually at least 7 degrees), making the approximation more robust. Meanwhile, another restriction about the intercept is used to control the range of gradients in early iterations to avoid early gradient exploding.

However, in the output layer, the activation function mainly depends on the type of prediction tasks. For instance, if the task is binary classification, then the sigmoid function will usually be selected, and in our experiment, even a 15-degree least squares fitting polynomial may fail to converge. Therefore, we propose a secure piecewise linear approximation to maintain the boundedness of the original activation function without the leakage of actual predictive values, whose function and pipeline are shown below:

${\displaystyle {\displaystyle y=\{\begin{array}{ll}{a}_{1}x+b_1& x<d_1\\ a_{i}x+b_i& x\in [d_{i-1}{d}_i],i\in [2,n-1]\\ a_{n}x+b_n& x\ge d_{n-1}\end{array}}}$(4)

First, inner products between the input values of the output layer and the output layer’s model parameters are computed under HE to obtain the inputs of activation function c(x), where one value may appear in multiple ciphertext slots according to the number of linear regions. Next, these calculated inputs are subtracted with discontinuous points ${d=d}_0,d_1,\dots ,d_{n-1}$ (d₀=0, and the order is disturbed randomly) and multiplied by large uniformly distributed noises ${e=e}_0,e_1,\dots ,e_{n-1}\mathrm{}$ sampled from $U(-a_i,a_i)$ (a_i is randomly selected from different distributions for each e_i), which are stored locally and have the same order of discontinuous points. Then, the resulting ciphertexts are decrypted and refreshed, where the refreshed ciphertexts $c(xe-de)$ and masked plaintexts $xe-de$ are obtained by the researcher. After obtaining masked plaintexts, the researcher generates interval indicator plaintexts ε with only 0 and 1 (here, 1 and 0 indicate positive and negative values in masked plaintexts, respectively), encrypts them (denoted by $c(\epsilon (xe-de))$), and sends these ciphertexts together with refreshed ones to corresponding data providers. Then, data providers also generate their binary indicator plaintexts $\epsilon (e)$ according to multiplicative noises and finally compute the piecewise linear approximation secretly as follows, which uses only two multiplication depths and can be easily implemented by “rotate and sum” operations:

${\displaystyle {\displaystyle \begin{array}{l}c(y)={\displaystyle \frac{a_1}{e_0}}c(x{e}_0)+b_1+{\displaystyle \sum _{i=1}^{n-1}{\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_i\cdot (2ϵ(e_i)-1)\cdot }\\ c\left(ϵ(x{e}_i-d_i{e}_i\right)+ϵ(e_i)-1)\cdot c(x{e}_i-d_i{e}_i)\end{array}}}$(5)

${\displaystyle {\displaystyle \begin{array}{c}slop{e}_i=\frac{a_{i+1}-a_i}{e_i}+b_{i+1}-b_i+d_i\left(a_{i+1}-a_i\right)\end{array}}}$(6)

‎

During the process, the only exposed information is masked inputs of the output layer’s activation function, from which an attacker is unlikely to infer any valid private data due to the zero-centered noises e and random order of $e(x-d)$. The researcher can also benefit from them by judging whether the absolute values of masked inputs have grown too large, which usually means gradient explosion, thus stopping the failed training in the early stage and significantly avoiding the waste of computational resources. The main limitations of this HE-based piecewise linear approximation are (1) the need for an extra round of communication, which can coincide with ciphertext refresh by designing the forward propagation algorithm of secure sliced-RNN (see more details in Multimedia Appendix 1 in the section “Secure RNN Example”) and (2) the low utilization rate of ciphertext slots, which is inversely proportional to the number of linear regions, but in many real-world medical deep learning analysis scenarios, it is still acceptable since the number of outputs in the output layer is usually much smaller than that in the hidden layer.

As for network structure, we apply a sliced-RNN architecture when the time step is large, which splits the long input sequence into multiple equal-length subsequences and uses a traditional RNN-like structure on each subsequence. Then, all hidden states at the end of each subsequence are treated as a long sequence, which is again equally partitioned with the RNN-like structure. The split-forward cycle will repeat several times until the final sequence length reaches the requirement for model prediction [42]. Due to its parallel short subsequences and multiple hidden layers, it has some advantages that may be more prominent in the HE environment. (1) Compared to traditional RNNs, because one hidden layer’s all subsequences can be calculated simultaneously, sliced networks are more parallelized (especially in lower layers), meaning that they can make use of ciphertext slots more efficiently in most cases and are better adapted to various kinds of parallelization acceleration as well as the newly proposed ciphertext matrix multiplication protocol. (2) Sliced-RNN typically can obtain more information about long-term relationships through multiple hidden layers with little computational cost, which is very important in computationally expensive HE-based scenarios. However, more complex network structures and more model parameters can also lead to overfitting problems in noisy data. Meanwhile, low-degree polynomials typically provide limited ability to approximate nonlinear relationships compared to nonpolynomial activation functions like ReLU and tanh, possibly causing slower convergence and underfitting problems. Therefore, the combination of them may have a good trade-off.

Results

Experimental Setup

Here, we use three main aspects to evaluate the performance of our proposed multiparty longitudinal patient data modeling protocol: (1) predictive accuracy: comparing accuracy among secure, centralized nonsecure, and single-site nonsecure RNN model training, including original and sliced network structure; (2) efficiency and scalability: the time needed to perform secure sliced-RNN model training and its growth with data size and training settings; (3) security: whether sensitive data and trained models are at risk of exposure.

Four datasets are used for the experiments, which are described below. The first comes from PhysioNet Challenge 2012 [44], containing records of 12,000 intensive care unit (ICU) stays of adult patients. Given the test values within the first 48 hours, the target is to predict in-hospital mortality. The second comes from PhysioNet Challenge 2019 [45], where we randomly select 8000 ICU stays with at least 48 rows of measurements. For every ICU stay, we only use the last 48 rows, and the target is to predict the occurrence of sepsis. The third and last one is obtained from the MIMIC-III dataset, where the former contains test values from 18,415 patients during the first 8 days of hospitalization and aims to predict the 28-day mortality, and the latter includes 48 minutes long densely collected information about blood pressure, heart rate, and respiration from 5465 patients and is used for prediction of acute hypotension episodes [43]. Each time step is 1 hour in the first two datasets, 1 day in the third dataset, and 1 minute in the last one, making them have 48, 48, 8, and 48 time steps, respectively. Finally, data standardization is performed on all datasets to transform the range of different features into the same scale, and for missing values, the last observation carried forward is applied if a previously observed value exists. Otherwise, mean imputation is used.

For the encryption parameter setting of approximated HE, the selection is shown below to achieve 128-bits of IND-CPA^D security and (123, 32)-bits of 2¹⁵-IND-CPA^D security, which match the most recent work about IND-CPA^D and the security standards proposed by HomomorphicEncryption.org group [40,55]. The degree of polynomial modulus is 8192, while the coefficient modulus is a product of 7 distinct 31-bit primes. The multiplication depth of the HE scheme is 4, so we choose a 31-bit-long integer as the scale factor. The key, error, and decryption distributions follow U{–1, 0, 1}, N(0, 3.2²), and N(0, 2^64.58), respectively.

For a real-world scenario, we implement our framework on different machines with Windows 10, each representing a data provider or researcher. The machines have AMD Ryzen 5 5600X CPUs running at 3.7 GHz, with 12 threads on 6 cores and 32 GB of RAM. Communication between any two machines is done through the Transmission Control Protocol with a bandwidth of 1 Gbps. To achieve the threshold variants of the CKKS scheme and support privacy-preserving multiparty sliced-RNN model training, we use the Microsoft Simple Encrypted Arithmetic Library v4.1.1 and have made some modifications [56].

Accuracy Comparison

We first evaluate the impact of ciphertext noises on the predictive accuracy of the trained model, where the noises are produced during encryption, decryption, and computational operations. Here, we use the PhysioNet Challenge 2012 dataset and a stratified random sample of 2/3 for the training set and 1/3 for the testing set. Given the same initialized model parameters and sliced-RNN configuration (number of features f=39, number of hidden layer and output layer’s dimensions $h,o=31$, mini-batch size $b=250$, learning rate $\alpha =0.001\times {0.5}^{0.01\times \mathrm{e}\mathrm{p}\mathrm{o}\mathrm{c}\mathrm{h}}$, time step setting $T=\left[\mathrm{48,2},\mathrm{2,2},\mathrm{2,3}\right]$, where the first number indicates the total time steps and the following ones mean the minimum subsequence of each hidden layer, number of epochs $ep=25$), the noise estimation and area under the curve (AUC) comparison are shown in Table 1. Here, both the inherent noise of ciphertexts and smudging noise have nearly no impact on model performance. Therefore, to accelerate the accuracy evaluation of our framework, we perform it in a pseudo-HE manner, which replaces all activation functions with approximated polynomial (ReLU, 2-degree for the first three datasets and 4-degree for the remaining one) and piecewise linear (sigmoid) functions, and still calculates in plaintext but limits all computation operations to addition, multiplication, and vector rotations.

Then, to demonstrate the model prediction performance of our proposed framework, we conduct a comparative experiment using all above-mentioned longitudinal patient datasets with (1) nonsecure RNN with ReLU or tanh in the hidden layer and all data centralized in 1 machine, where ReLU is closer in shape to our low-degree polynomials, and tanh is more commonly used in RNN, (2) nonsecure sliced-RNN with ReLU in hidden layers and all data centralized in 1 machine, and (3) secure RNN with data distributed in multiple machines with replaced polynomial activation functions like ours. The training settings are shown in Table 2, which repeats 50 times using 10 different random seeds for the stratified sampling of the training and testing set, and for each stratified sampling, 5 different random seeds are used for model parameter initialization (all use the same random seeds in one experiment). Here, to achieve faster convergence and avoid the gradient explosion problem as much as possible, we choose Xavier initialization for parameters of the hidden layer in RNN (typically, He initialization should be used on the ReLU function, but in our experiment, it is significantly slower to converge and more vulnerable to gradient explosion in RNN) and output layer, and He initialization for parameters of the hidden layer in sliced-RNN. Class weights are also used to address the label imbalance of all datasets. Finally, after each epoch, the overall loss of the testing set is calculated secretly and decrypted for possible early stops, which usually obtain AUCs slightly lower than the theoretical maximum AUC.

The experimental results are shown in Figure 5 and Table 3, where the former shows the variation of maximum AUC with the number of epochs for the proposed and compared solutions, and the latter lists their prediction performance and the corresponding number of epochs required to reach the minimum overall loss in the testing set (averaged over 50 repetitions of the experiments with different random seeds). In general, the use of low-degree approximated functions does lead to more epochs for near-optimal prediction performance, but surprisingly, the model predictive ability of secure solutions is overall a little better than nonsecure solutions, where the performance is relatively much closer in between nonsecure and secure RNNs and significantly different between nonsecure sliced-RNN and ours. The uncommonly low maximum AUC of the nonsecure sliced-RNN model may be related to the fact that its model structure is overly complex compared to traditional RNNs, while longitudinal clinical data tend to have much more missing and biased data than many other fields due to the unsuccessful measurements and specificity of patients during long observation windows (especially in the first 3 datasets since the interval between time steps is much longer than that in the last dataset) [57], which together lead to overlearning of the training set’s noisy features and severe overfitting. By contrast, our solution significantly mitigates the overlearning problem of the training set by limiting the ability to represent nonlinear relationships through low-degree polynomial activation functions, making it more robust against noisy features, thus obtaining even higher prediction performance on all 4 datasets compared to the original nonsecure RNN. It is also the first case, to the best of our knowledge, where the HE-based trained model has better predictive ability than the corresponding nonsecure one with the same training data, while similar previous studies have shown either significant or minor losses in accuracy [30,36].

Compared to the secure RNN, our solution also provides better training results overall, significantly outperforming it in the 1st, 2nd, and 4th datasets, either requiring only much fewer epochs while having similar prediction performance (1st dataset) or having explicitly better prediction performance with relatively closer epochs (2nd and 4th datasets). Notably, the 4th dataset has much shorter testing time intervals and fewer missing values than the other 3 datasets, in which secure RNN obtains the worst prediction performance, showing its limited ability to handle nonlinear relationships with less noisy features. In contrast, ours achieves a better trade-off and obtains the best prediction performance. On the other hand, the secure RNN’s trained model performs slightly better in the third dataset, which may be caused by the dataset’s much fewer time steps. Such short sequences are insufficient to fully utilize the long-term relationship acquisition advantage of sliced-RNN, while the sliced structure impairs neighboring relationships to some extent between parts of time steps.

Furthermore, to demonstrate the advantage of cross-institutional collaborative analysis over single-site analysis, we conduct a comparative performance evaluation between our secure sliced-RNN models (trained on the entire dataset) and nonsecure RNN models (trained on only a partial dataset) with the same test sets in Table 3. We simulated a realistic scenario where data are equally partitioned between two institutions (A and B) but with an imbalanced label distribution (positive rate ratio of 2:3). As shown in Table 4, our secure sliced-RNN model, trained on the entire dataset, demonstrates significantly superior predictive performance across all 4 medical datasets compared to models trained on isolated institutional data. Critically, this advantage is consistently observed under both label imbalance and balanced conditions. In the first 3 datasets, where positive cases represent only 10%-20% of the population, our model’s performance gain is occurred most in the AUC and recall—precisely the metrics most crucial for identifying rare but clinically critical events. For instance, on the PhysioNet Challenge 2012 dataset, SRNN (sliced recurrent neural network)-secure achieves a recall of 0.7981, dramatically outperforming both RNN A (0.5937) and RNN B (0.6057), demonstrating a markedly enhanced ability to identify true positive cases without being hindered by the class imbalance. This superior performance extends to the fourth dataset, which exhibits a balanced label condition. Here, SRNN-secure not only achieved the highest AUC but also excelled across all metrics including accuracy, recall, and F₁-score, showcasing its robustness and ability to integrate diverse data patterns effectively.

‎

Notably, these results are achieved under a constrained simulation with only two institutions. We anticipate that this performance gap would be substantially wider in real-world multi-institutional collaborations involving more data providers with greater heterogeneity. Our framework effectively harnesses such diversity to build more robust and generalizable models, which is unattainable by any single institution working in isolation.

Analysis of Long-Term Training Stability

The use of low-degree polynomial approximations (eg, quadratic functions for ReLU) raises concerns about gradient explosion over extended time steps, as the hidden input may exceed the fitting interval during the long-term forward and backward propagation, making the approximations completely unreliable. Therefore, to empirically validate stability, we conducted a stress test on the PhysioNet Challenge 2012 dataset by artificially extending the sequence length from 48 to 96 and 192 steps via linear interpolation of the original 48-step trajectories. Across 50 independent runs for each sequence length (the same random settings of Table 2, shown in Multimedia Appendix 1 in the section “Long-Term Training Stability With Approximated Activations”), we find that (1) no hidden input exceeds the fitting interval, even at 192 steps; (2) both the hidden input range and the $L_{\infty }$-norm of gradients show no significant increasing trend as sequence length grows (48-step: hidden input range [−41.560, 39.474], $L_{\infty }$-norm 0.0147; 192-step: hidden input range [−43.692, 43.297], $L_{\infty }$-norm 0.0142). These results confirm that our quadratic ReLU approximation provides stable training for sequence lengths up to at least 4 times the original 48-step limit. We attribute this stability to (1) the sliced-RNN architecture, which effectively limits the propagation depth; (2) the bounded nature of the input features after standardization; and (3) the smoothing effect of low-degree polynomial approximation, which reduces sensitivity to noise and mitigates overfitting.

Time Consumption and Scalability

The scalability of our protocol is shown in Figure 6, where we test our framework in different experimental settings. Specifically, as the number of time steps grows, both the computation and communication time increase proportionally; as h and o grow, the computation and communication time increase near quadratically and linearly, respectively, since matrix multiplication and outer product computation form the majority of computational overhead and most weight matrices enlarge quadratically with hidden and output dimensions; in contrast, as the number of data providers grows (each of them maintains the same local mini-batch size), the communication time increases linearly, but each data provider’s computation time varies very little. These results are consistent with previous similar HE-based works [12,36], meaning that adopting the sliced structure does not negatively impact the scalability of time-series data model training. Besides, the lightweight encryption parameters and low multiplication depth of the fitting polynomials and matrix multiplication in our framework also effectively reduce computational overhead, with about 1 minute to perform one iteration of model training in real-world longitudinal patient data. In comparison, in Sav’s work [36], 7-degree polynomial fitting tanh and ciphertext matrix multiplication with 3 multiplication depths are used, which leads to much larger encryption parameters and operation times of ciphertext multiplication, thus requiring nearly the same time to perform one iteration of RNN model training in only much smaller settings in their experimental results (the same hidden dimensions, doubled local mini-batch size, and 4 time steps).

‎

Furthermore, we conducted a simulation-based study to emulate various wide area network (WAN) conditions and estimate the impact of communication latency in our framework. Specifically, we artificially imposed three levels of round-trip time (RTT: 20, 50, and 100 ms) and evaluated the additional training time under different numbers of time steps and data providers for each latency setting. The results are presented in Multimedia Appendix 1 in the section “Impact on Communication Latency.” Notably, the additional training time increases linearly with the latency, and the growth is much less significant with the increasing of time steps and data providers. This is because the number of ciphertext communication rounds scales only logarithmically with the number of time steps due to the sliced network structure, and remains constant with the number of data providers, as the ciphertext refresh and partial decryption protocol requires a fixed number of communication rounds independent of party count. These findings demonstrate that our framework maintains good scalability under real-world WAN conditions and remains practical even under relatively high latency (eg, 100 ms), as the additional training time constitutes only a small fraction (approximately 10% in the default setting) of the total training time under nearly zero-latency conditions.

Security Analysis

The security of our proposed framework is discussed in the following 4 aspects.

Secret Keys Held by Each Party

We consider two situations for using secret keys. One is collective key generation, and the other is partial and final decryption, where the security of secret keys in the former case has already been proved by the simulation paradigm with all but one colluding party [58]. In the latter case, the risk of secret keys arises mainly from (1) the possible exposure of secret keys due to the delivery of partial decryption shares and (2) a full key recovery attack specialized for approximate RLWE-based HE that is mainly initiated by the final decrypting party [40], both of which can be solved by adding large smudging noises to achieve a sufficient level of statistical security. In addition, when the number of ciphertext refreshes and decryption reaches a certain value, all keys will be regenerated to limit the decryption query times under every set of keys, thus ensuring IND-CPA^D.

The noise estimation is non-worst-case, while Guo’s [37] attack showed that, only with worst-case noise estimation, the noise-flooding countermeasures can achieve IND-CPA^D. However, in their adversary setting, the attacker can select arbitrary inputs and HE computation protocols, which is not allowed in passive adversaries [59,60]. Therefore, this attack is not applicable to our framework, and non-worst-case noise estimation is sufficient to overwhelm ciphertext noise and plaintext information (54.77 bits vs 22.48 and 35.64 bits, shown in Table 1).

Whether the Researcher Is Honest (Security of Model)

All data providers can only obtain masked intermediate results during the whole model training process, where the researcher partially generates the smudging noises. As long as the researcher’s noise shares are large enough (the size of noise added by each party during the ciphertext refresh process has nearly no impact on the accuracy of subsequent computation, so it is easy to achieve), an attacker will not be able to infer any valid training model information from them. Meanwhile, for initial model parameters, the initialization is performed by all parties, meaning that they cannot be directly exposed even with all corrupted data providers.

Whether a Data Provider Is Honest (Security of Its Data)

For an honest data provider, the only exposed information is noisy logit values, which are eventually decrypted by the researcher and used to determine the belonging intervals for piecewise linear fitting functions secretly. Since the zero-centered multiplicative noises and the order of discontinuous points are randomly generated only by the honest data provider, an attacker is unlikely to extract the exact logit values. Besides, the honest data provider can disrupt the order of data training for each epoch, making it very hard for an attacker to match multiple logit values of the same sample over multiple epochs, thus further avoiding the implied data leakage from such aggregated information.

Extension to Active Adversaries

While our current framework assumes passive adversaries, its design is compatible with extensions to handle active adversaries—who may deviate from the protocol by submitting manipulated gradients or false partial decryption shares. Recent advances have demonstrated the feasibility of combining fully HE with Byzantine-robust aggregation to defend against such malicious behavior. For example, Lancelot introduces a mask-based encrypted sorting mechanism that enables robust aggregation rules such as Krum and median to operate directly on encrypted model updates, effectively protecting against malicious clients that attempt to poison the global model [61]. It achieves this enhanced security with little computational overhead, taking only 18× the time of plaintext operations (compared to 384× for the open-source fully homomorphic encryption [OpenFHE] approach). Therefore, it can be integrated into our framework to provide cryptographic privacy and Byzantine resilience simultaneously without incurring prohibitive computational or communication overhead.

Discussion

Principal Findings

The above experimental results demonstrate that our proposed framework successfully achieves a practical trade-off among security, accuracy, and efficiency. First, it enables secure multiparty RNN training on distributed longitudinal patient data without exposing raw data or model parameters. Second, it maintains good stability in long-term datasets, shows significant accuracy improvement compared to simulated single-site analysis, and achieves predictive performance comparable to—in some cases even exceeding—nonsecure centralized models (eg, AUC of 0.8480 vs 0.8404 on the PhysioNet Challenge 2012 dataset). Third, it maintains practical efficiency, requiring approximately one minute per training iteration under realistic settings, with favorable scalability as data size and number of participants increase. These findings confirm the feasibility of privacy-preserving deep learning on medical time-series data using threshold HE.

Interpretation, Implications, and Comparison With Prior Work

The integration of threshold approximated HE with noise engineering shows excellent potential. The large noises added by various parties during ciphertext refresh and decryption, together with the underlying plaintext message, effectively form a linear secret sharing scheme such as Mouchet’s collective bootstrapping protocol in the threshold BFV scheme, ensuring the security of the exposed masked values [58,62]. More importantly, we leverage this secret-shared plaintext domain to perform specific linear operations—namely, data scaling and encoding changes—that would otherwise require expensive homomorphic operations (specifically, ciphertext-plaintext multiplication and ciphertext rotation). This strategic optimization further enhances the multiplicative depth availability and overall ciphertext computation efficiency. Consequently, under a computable ciphertext environment, substantial noise can be injected into intermediate and final results during refresh and decryption without compromising computational precision. This approach guarantees the complete nonrecoverability of any exposed information and provides robust defense against statistical attacks targeting both threshold and approximated HE. In contrast, when relying solely on differential privacy, a challenging trade-off between security and accuracy invariably persists, particularly within RNN-like structures [26].

Our comparison of predictive ability between secure and nonsecure solutions also shows the promise of the proposed framework for high-quality multiparty privacy-preserving analysis on datasets with noisy and large parts of missing features, which is common in the biomedical field. Almost all previous HE-based machine learning and deep learning model training works have focused on fitting polynomials that are as close as possible to the original activation functions, mainly validated on high-integrity publicly available datasets, and expected to obtain analytical results that are either nearly the same as the corresponding nonsecure central solution or having an acceptable loss of accuracy with the same or a little larger number of iterations [12,30,36]. In contrast, we use more complicated networks and low-degree functions to fit the ReLU function with a large fitting interval, introducing a larger fitting noise, but achieving higher AUCs instead. This finding indicates that future HE-based iterative computing research may not only be simple mimics of existing algorithms but rather a comprehensive optimization problem combining data, polynomial activation function selections, network structure, and neuron structure, which can be very challenging and needs further exploration.

The significant difference between initial settings of maximum epoch times and final experimental obtained early stop points also indicates the importance of exposing specific information in HE-based iterative computing solutions due to the time-consuming and invisible nature of the calculation. For example, after a certain number of epochs, a privacy-preserving inference should be made on the testing set (well suited to our proposed matrix multiplication, as the number of testing samples is usually much larger than the batch size), and then the plaintext of some testing set’s aggregated information should be obtained by the researcher to facilitate timely stop of training, rather than empirically specifying hyperparameters, and suffering useless model training, or even overfitted results, leading to waste of a large amount of computational, storage, and network resources. Here, for a better trade-off among efficiency, security, and utility, we select overall loss as the indicator since it is straightforward to compute under an HE environment, very close to the theoretical best prediction performance (at least in our secure sliced-RNN, as shown in Figure 5), and reveals only one value per computation, which is hard for attackers to infer private information. On the contrary, if we directly use prediction performance (like AUC) as an indicator, the computation in the HE environment will be much more complex and impractical because of the limitations of ciphertext operations and efficiency.

Another advantage of our framework is the potential to adapt well to longitudinal patient data with different overfitting risks and task types without additional overfitting prevention settings. In comparison, HE-based solutions with perfectly fitted activation functions (as nonsecure ones, for example, using a proposed piecewise linear function to perform ReLU) are more or less affected by overfitting and can only obtain comparable prediction performance to ours when using overfitting prevention measures (such as L2 regularization) with proper hyperparameters at the cost of obviously more epoch times overall (see more details in Multimedia Appendix 1 in the section “Performance With L2 Regularization”). This is essential in privacy-preserving real-world medical data analysis scenarios, where researchers often do not have direct access to research data and, therefore, have difficulty determining overfitting prevention hyperparameters. Federated tuning methods have problems with privacy and heterogeneity: private information leakage mainly comes from sharing of private local model parameters, weights, and optimal hyperparameters [20,21,63], while heterogeneity makes it difficult to approach the real global optimum and is quite common in the medical research domain since different medical institutions often have different treatment abilities, scopes of patients, and data collection strategies, causing distinct feature, label, quantity, and noise distributions among all data providers [64,65]. Centralized tuning methods like grid search and Bayesian optimization are accurate, but the operation time may grow very large with extra hyperparameters due to their iterative nature and the expensive computational cost of HE operations. Therefore, with fewer hyperparameters, our solution may save substantial time when used in the real world.

Limitations

Here, we test our framework in many-to-one scenarios, but in the real world there are quite a few many-to-many scenarios as well (eg, prediction of acute exacerbation at fixed time intervals during hospitalization), where if we continue using the proposed secure piecewise linear approximation, an excessive number of ciphertexts may be produced to lower the computational efficiency due to the low utilization of ciphertext slots. In this case, some alternative polynomial fitting techniques exist for large intervals. For instance, the domain extension algorithm by Cheon et al [66] can extend the valid fitting interval of an approximated polynomial by a factor of 2.45 for every two multiplicative depths consumed, which is very useful when the original polynomial is high-degree (eg, for sigmoid, 9-degree or higher), thus requiring larger encryption parameters or more rounds of ciphertext refresh. As for ciphertext matrix multiplication, Jang et al [33] recently proposed a variant of the CKKS scheme called MatHEAAN, which directly supports ciphertext rotations in a matrix manner and occupies only two multiplication depths without losing any ciphertext slots. However, as it is not yet integrated with the RNS, its current efficiency does not surpass our RNS-based scheme, though future advancements may make it a viable alternative. Besides, while our stability analysis demonstrates robustness up to 192 steps, we acknowledge that the interpolated sequences may not fully capture the statistical properties of real clinical time series of such lengths, which will be further evaluated in the future.

Conclusions

This paper proposes a protocol for privacy-preserving multiparty RNN analysis of longitudinal patient data using threshold HE. The experimental results demonstrate that secure training can achieve even higher predictive accuracy compared to nonsecure centralized models in some cases, with practical efficiency and good scalability. More broadly, these findings challenge the prevailing assumption that cryptographic privacy inevitably sacrifices accuracy, and instead reveal an unexpected synergy: low-degree polynomial approximations, chosen primarily for computational efficiency, also act as an implicit regularizer that mitigates overfitting on noisy clinical data. This insight points toward a new design paradigm for privacy-preserving medical AI—one where cryptographic protocols and model architectures are co-designed rather than naively combined. We believe that with existing clinical information integration platforms, our framework will promote multiparty longitudinal medical research and help accelerate knowledge discoveries.

Some future works still exist to make it more practical. First, GPU (graphics processing unit) acceleration can significantly improve the efficiency of basic ciphertext operations, while the software pipelining technique can better accelerate our framework than traditional RNNs since different subsequences can be executed in parallel under a sliced structure. Second, consider the use of more complex neuron structures and apply them to research scenarios where more data providers provide the longitudinal patient data, then explore the predictive performance of ciphertext computation solutions in different configurations on data with a varying degree of missing and noisy features. Furthermore, it is also essential to evaluate the impact of hyperparameter tuning with different heterogeneity in real-world medical data distribution scenes. If needed, an efficient privacy-preserving tuning method may be further researched.